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Rapid inversion of angular deflection data for 
certain axisymmetric refractive index distributions 


R. Rubinstein and P. S. Greenberg 


Certain functions useful for representing axisymmetric refractive-index distributions are shown to have 
exact solutions for Abel transformation of the resulting angular deflection data. 
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The tomographic reconstruction of axisymmetric re- 
fractive-index distributions from angular deflection 
data (e.g., schlieren, moire) is a subject of consider- 
able interest, particularly in the areas of combustion 
science and fluid physics. Although the mathemati- 
cal solution to this problem was presented by Abel 1 
more than a century ago, the formal solution is 
seldom used in practice, primarily because of the 
singular behavior of the integrand at r = 0 and the 
requirement for spatially differentiated, often noisy, 
experimental data. The latter is less of a concern in 
the present context, since the resulting angular devia- 
tions are physically associated with the spatial gradi- 
ents of the field itself; thus one avoids problems 
associated with the differentiation of noisy data. A 
number of schemes exist for the numerical inversion 
of line-of-sight data, and their properties have been 
extensively investigated and compared. 2 We show 
here that certain classes of function suitable for the 
representation of axisymmetric refractive-index fields 
have exact solutions with respect to the Abel inver- 
sion of the resulting angular deflection data. These 
solutions include both Gaussian and Hermite polyno- 
mials, the latter being of particular value since Her- 
mite polynomials form an orthonormal basis for 
radial decomposition. The solutions are physically 
useful owing to their asymptotic behavior as r — > oo. 
Since the solutions are in the form of combinations of 
known functions and recurrence relations, the inver- 
sion reduces to a trivial lookup table operation involv- 
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ing precomputed coefficients. This process can be 
implemented extremely rapidly, even when high poly- 
nomial orders are retained. In addition the process 
of basis decomposition with a finite cutoff effectively 
smooths the data, which diminishes the effects of 
spatial noise. 

When a planar refractive-index distribution 5 is 
axisymmetric, 5 = 8(r), and satisfies the conditions 

f l d5 l 

J |dr| ds ^ 1] ( la > 

8(r) •-» 0 for r -> oo, (lb) 

where L denotes any ray path and ds is the correspond- 
ing element of arc length, then the optical path length 
and the angular deflection along any ray can be 
evaluated as if the ray did not depart from its initially 
linear trajectory. Here we use the standard conven- 
tion 8 — n - 1, where n is the total value of the 
refractive index. Condition (lb) means that we de- 
note by 8 the departure of the refractive index from 
its value in empty space. This approximation leads 
to the well-known relationships 3 between angular 
deflection (e.g., schlieren) data e(y) and the refractive- 
index distribution 8(r): 

* d5 dr 

dr (r 2 - y 2 ) 1/2 ’ ^ 2a) 

8(r)=-i| e(y) (y2 J^. 2)1/2 • (2b) 

Equation (2a) represents the direct problem of evalu- 
ating the angular deflection, and Eq. (2b) the inverse 
problem. Figure 1 illustrates the coordinate system 
of the initial ray trajectories and their resulting 
deflections. 
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Fig. 1. Geometry of ray deflections in an axisymmetric refractive- 
index distribution n(r): y, entrance coordinates of the initial rays; 
ds, differential length along any given ray; r, radial coordinate of 
the refractive-index distribution; e(y 2 )> exit angle of the ray originat- 
ing at the coordinate y l = y 2 * 


A simple pair of functions related by Eqs. (2) is 
J[ v exp(-v 2 ) — 2 _ d ^ - 2 — 2 = l (tt) 1 2 exp (-r 2 ), (3) 

which follows from elementary changes of the vari- 
able in the integral for the gamma function. This 
equation determines the angular data corresponding 
to a Gaussian refractive-index distribution, 

8(r) = A exp (— r 2 ), (4) 

or the refractive-index distribution when the angular 
data are of the form 


e(y) = By exp(-y 2 ). (5) 

Substitution of either Eq. (4) or (5) into Eq. (2a) or 
(2b), respectively, reduces at once to Eq. (3). 

These results can be generalized by 


f 


v 2n+1 exp(-v 2 ) 


dv 

(v 2 - X 2 ) 1 '* 


= o (*n , ) 1/2 JC 2n+1 exp( —x 2 


n + ~>x 2 


( 6 ) 


where U is a solution of Rummer’s equation. 4 We 
can evaluate the function U for n , an integer, by 
applying the recurrence relations 



and the initial value 



(tt) 1 2 exp(z 2 )erfc(z). 


( 8 ) 


From Eqs. (6) and (7) follow 



These results generalize the stated results above as 
follows: Define a Gaussian polynomial to be the 

product of a polynomial and the Gaussian function 
exp(-r 2 ). If the refractive-index distribution is an 
even Gaussian polynomial, 

8(r) = (a 0 4* a 2 r 2 + * * -)exp(-r 2 ), (10) 

its derivative is an odd Gaussian polynomial, and the 
corresponding angular data are the odd Gaussian 
polynomial found by substitution of the results of 
Eqs. (6)-(9) into Eq. (2a). Conversely, if the angular 
data are a given odd Gaussian polynomial, the corre- 
sponding even Gaussian polynomial refractive-index 
distribution follows by substitution of the results of 
Eqs. (6)-(9) in Eq. (2b). 

These exact results imply a simple approximation 
scheme when the angular data can be accurately 
represented by an odd Gaussian polynomial; because 
of the asymptotic return of the refractive index to its 
reference value in many combustion science and fluid 
physics applications, this approximation is reasonable. 
It is necessary first to approximate the angular data 
by an odd Gaussian polynomial. The orthogonality 
of Hermite polynomials makes this process straight- 
forward, namely, 

z(y) ~ 2 aiHj{y)exp{-y 2 ) ( 11 ) 

i odd 


with 


= r w 

J_» I-*) 1 2 i!2‘ ‘ 

Inspection of Eqs. (2) demonstrates that the refractive- 
index distribution is given by an even Gaussian 
polynomial. 

The recurrence relation [Eqs. (7)] can be applied to 
expedite the calculation. Since the right-hand side 
of Eq. (6) is a Gaussian polynomial, we define 

= vv,* 2 ). 

Then Eqs. (7) provide the recurrence 

V(n, x 2 ) = [x 2 + n - ~ 1- x 2 ) 

+ x 2 (l - n)V(n - 2, x 2 ), (12) 



We define H n m by 

H„(x) = 2 H n , m x m . 

m 2:0 
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Explicitly 

H 0 0 = 1, H 0 m = 0 for m > 1, 

H lt0 = 0, Hu = 2, H l m = 0 for m > 2, 
tf n ,o = -2 (n - 1)H„_ 2 ,o for n > 2, 

= “2(n ~ 

for n > 2, m > 1. (13) 

Then, writing 

e(y) = 2 a n H n (y)exp(-y 2 ) 

n odd 

= X o n f/ n>+ iy^ +1 exp(-y 2 ), 

n odd,p^0 

we give the corresponding refractive index as 


n (r)=(l/ 2 ){e-<r-i )«-2 - e -(r+ 0 ** 2 } 
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Fig. 2. Error in seven-term reconstruction of a model refractive- 
index field. 


5(r) = 


1 p dy 

w J r ( y 2 ~r 2 ) 


2 a n H nt2p+1 y^ +l exp(-y 2 ) 

n odd,p>0 


i |2p + 1 


2(TT) 1 2 n odd,p > 0 


with C independent of n. Therefore noise in the in- 
version is then bounded by noise in the data, specifi- 
cally the magnitude of the noise does not grow if the 
number of terms used in the inversion is increased. 

Note the equality 


i,2p+ 1 


x exp{-v 2 )ul^P + 2 ’ 

= -^ 2 a n H n , 2p+1 exp(-v 2 Mp > v 2 ) ( 


n odd,p ^0 


f 

Jo 


v 2n+! e xp( V 2 )dv = g n1 * 


and the estimate 


which we can rapidly and accurately evaluate to any 
order by using Eqs. (12) and (13). 

These results are stated here in terms of dimension- 
less variables. In inverting experimental data the 
variables should be scaled so that the decay of the 
Hermite polynomials appropriately matches the de- 
cay of the data. 

As a representative application we reconstructed 
the refractive-index distribution 


J v 2n exp( — v 2 )dv ^ n\. 
o 

Let 

V(n, x 2 ) = a 0 x 2n + a x x 2n ~ 2 + • • • + a 2n . 

Then the recurrence formula given by Eq. (12) pro- 
vides the bound 


8 (r) = exp[-(r - l) 2 ] + exp[-(r + l) 2 ], 

using seven-term Hermite interpolation of correspond- 
ing numerically evaluated angular data. The results 
in Fig. 2 show that the reconstruction is very accurate. 
One can reduce the error in the reconstruction fur- 
ther simply by taking higher-order Hermite approxi- 
mations of the data. Because the approximation 
requires functions that are evaluated by simple recur- 
re nce relations and coefficients that are evaluated by 
simple numerical quadratures, the computer re- 
sources required by the inversion are negligible. 

An advantage of this procedure over direct numeri- 
cal Abel inversion is that least-squares curve fitting is 
a smoothing process that reduces the noise sensitivity 
of the computation. 6 To bound the effect of noise on 
this inversion, we establish the inequality 


| V(n , x 2 )exp( -r 2 ) | dr < C v 2n+1 exp(-v 2 )dv 

Jo 


|a„| < 3(2p - 1)!!/2 P+1 ~ 3p! forp ;> 1, 

where as usual (2p — 1)!! = (2 p - l)(2p 3), . . . , 1. 

Therefore 


f 

Jo 


|V(n, x 2 )exp(-x 2 )dx j J 

< 6tt 1/2 ' 


p 2 n + l exp(“ v 2 )dv 


v (gp - 1 ) 1 -! ( 2 ” - ?p - 11 !! , 1 

,Ti 2P+1 


n\2 2n ~ 2p 


6TT 


. 1/2 


X p!(n -p)!/n! + 1 

pal 


The inequality 

g 

2 pK« _ p)!/ n - ^ 3 

0 <.p<>n 

then establishes inequality (14). Inequality (15) has 
been evaluated numerically and is shown to peak at 
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its maximum value for n = 3 and 4, from which point 
it decreases monotonically with increasing n. 
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